% @author:夜剑听雨
% E-mail:2530595378@qq.com
% 2025.08.29
% 反演主程序
clc;clear;close all;

%% 加载观测的频散曲线
load('./save_dir/pick/shot3_picked.mat');

Freq = (picked_data.frequency)';        % 获取频率向量 (Hz)
pvObserved = picked_data.velocity1;  % 拾取的基阶频散曲线
layNum = 4;    % 地层数量
% 纵波速度、地层密度，反演过程中保持不变
Vp = [1000 1000 1000 1000];       % 纵波速度(m/s)
Den = [2000 2000 2000 2000];      % 地层密度(kg/m3)

% 写入结构体ObsInfo
ObsInfo = struct(...
    'layNum', layNum, ...
    'freq', Freq, ...
    'vp', Vp, ...
    'den', Den, ...
    'pvObserved', pvObserved, ...
    'misfitType', 'L2');

%% Parameter
% 设置VS和层厚h的上界和下届
% 模型1上下边界
% Vs_ub = [300 400 500 800];h_ub = [3 5 8];
% Vs_lb = [100 200 300 500];h_lb = [1 2 4];

% 模型2上下边界
% Vs_ub = [300 300 500 600];h_ub = [6 3 6];
% Vs_lb = [100 100 300 400];h_lb = [2 1 2];
% 模型3上下边界
Vs_ub = [150 300 600 500];h_ub = [6 3 6];
Vs_lb = [50 100 400 300];h_lb = [2 1 2];

ub = [Vs_ub h_ub];
lb = [Vs_lb h_lb];
maxIter = 100;  % 反演最大迭代次数
numPop = 50;    % 反演中初始化变量的数目
% 写入结构体InvParameter
InvParameter = struct(...
    'ub', ub, ...
    'lb', lb, ...
    'maxIter', maxIter, ...
    'numPop', numPop ...
    );

tic % 反演开始计时
loopNum = 5;
Loss = zeros(maxIter+1 ,loopNum);
Ini_model = zeros(loopNum, 2*layNum-1);
Inv_model = zeros(loopNum, 2*layNum-1);
parfor k=1:loopNum
    disp("第" + num2str(k) + "次反演开始！")
    % 反演核心程序
    [loss_out,model_out]=PSOInv(ObsInfo, InvParameter);
    % 保存反演结果
    Loss(:,k) = loss_out;
    Ini_model(k,:) = model_out(1,:);
    Inv_model(k,:) = model_out(end,:);

    disp("第" + num2str(k) + "次反演结束！！！")
end
time1 = toc;
disp("反演结束！用时：" + time1)

% 保存反演数据
outfolder = './Save_dir/Inv'; % 指定想要创建的文件夹名称及其路径
if ~exist(outfolder, 'dir') % 检查文件夹是否存在
    mkdir(outfolder); % 如果文件夹不存在，则创建文件夹
end

save(fullfile(outfolder, 'PSO_loss3.mat'), 'Loss'); % 保存文件
save(fullfile(outfolder, 'PSO_Ini_model3.mat'), 'Ini_model'); % 保存文件
save(fullfile(outfolder, 'PSO_Inv_model3.mat'), 'Inv_model'); % 保存文件

